library(tidyverse)
library(lubridate)
library(viridis)
library(sf)
library(scales)
chi <- read_csv("raw-data/crime_data/Chicago/chicago_crimes_2018.csv")
## create date/time variables

# date ony variable
chi$Date_simp <- substr(chi$Date, 1, 10)
a<-as.factor(chi$Date_simp)
abis<-strptime(a,format="%m/%d/%Y", tz="UTC") #defining what is the original format of your date
chi$Date_simp<-as.Date(abis,format="%Y-%m-%d")

# time variable
chi$time <- format(strptime(trimws(substr(chi$Date, 11, 22), which="both"), "%I:%M:%S %p"), "%H:%M:%S")

Chicago Crime 2018 Descriptive Statistics

Top overall crimes

# top overall crimes
top_crimes <- chi %>%
  select(Date, `Primary Type`) %>%
  count(`Primary Type`) %>%
  arrange(desc(n)) 
top_crimes
## # A tibble: 32 x 2
##    `Primary Type`          n
##    <chr>               <int>
##  1 THEFT               64946
##  2 BATTERY             49772
##  3 CRIMINAL DAMAGE     27792
##  4 ASSAULT             20369
##  5 DECEPTIVE PRACTICE  18368
##  6 OTHER OFFENSE       17063
##  7 NARCOTICS           12841
##  8 BURGLARY            11708
##  9 MOTOR VEHICLE THEFT  9992
## 10 ROBBERY              9689
## # … with 22 more rows

Top “night” crimes

# top "night" crimes
top_night_crimes <- chi %>%
  select(Date_simp, time, `Primary Type`) %>%
  filter(time >="22:00:00" | time <= "06:00:00") %>%
  count(`Primary Type`) %>%
  arrange(desc(n)) 
top_night_crimes
## # A tibble: 31 x 2
##    `Primary Type`          n
##    <chr>               <int>
##  1 BATTERY             16249
##  2 THEFT               12582
##  3 CRIMINAL DAMAGE      9151
##  4 ASSAULT              4031
##  5 ROBBERY              3522
##  6 OTHER OFFENSE        3360
##  7 BURGLARY             3064
##  8 DECEPTIVE PRACTICE   2940
##  9 MOTOR VEHICLE THEFT  2832
## 10 WEAPONS VIOLATION    1885
## # … with 21 more rows
## plot primary crime types of time
top_five<- top_crimes %>%
  head(5) %>% pull(`Primary Type`)

# chi %>%
#   select(Date_simp, `Primary Type`) %>%
#   filter(`Primary Type` %in% c(top_five, "NARCOTICS")) %>%
#   arrange(Date_simp) %>%
#   group_by(Date_simp) %>%
#   count(`Primary Type`) %>%
#   ggplot(., aes(Date_simp, n, color=`Primary Type`)) +
#   geom_line() + theme_classic()
#### theft crime counts overtime, 2018
# date variables for geom_tile plot
chi$month <- month(as.POSIXlt(chi$Date_simp, format="%Y/%m/%d"))
chi$year <- year(as.POSIXlt(chi$Date_simp, format="%Y/%m/%d"))
chi <- chi %>%
  mutate(month_year = paste0(as.character(month),"-",as.character(year)))
chi$weekday <- wday(
  chi$Date_simp,
  label = T,
  abbr = TRUE,
  week_start = getOption("lubridate.week.start", 1)
)
chi$weeknum <- week(chi$Date_simp)
chi$daynum <-day(chi$Date_simp)
chi %>%
  # select(Date_simp, `Primary Type`) %>%
  filter(`Primary Type`=="THEFT") %>%
  arrange(Date_simp) %>%
  group_by(Date_simp) %>%
  mutate(n = n()) %>%
  ggplot(., aes(x = weeknum, y = weekday, fill = n)) +
  geom_tile(col = "black", width = .95, height = .9) +
  # facet_wrap(vars(year), strip.position = "top", nrow = 3) +
  scale_fill_viridis(option="inferno", #end = .9, labels = function(x) paste0(x, "%"),
                     name="", direction = -1) +
  scale_x_continuous(expand = c(0, 0),
                    breaks = seq(1, 52, length = 12),
                    labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
                               "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))  +
  xlab("") + ylab("") +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        text = element_text(family="Lato"),
        legend.justification = c(.95, 0),
        legend.position = "bottom",
        legend.margin = margin(-0.25,-0.5,0,-0.5, unit="cm"),
        plot.margin =  margin(.5,.5,.5,.5, unit="cm"),
        legend.key.width = unit(.8, "cm"),
        plot.title = element_text(vjust=2.4, size=14),
        plot.subtitle = element_text(size=9, vjust=1),
        plot.caption = element_text(size=7, vjust=-15)
        ) +
  ggtitle("Chicago Theft Crime Counts, 2018")

# load community area shapefile
commun_area <- read_sf("raw-data/Boundaries-Community-Areas_current/geo_export_14339f9c-fd8b-4301-8af0-035b3fbf2c00.shp")
## https://datahub.cmap.illinois.gov/dataset/community-data-snapshots-raw-data/resource/96bc2e7d-9276-4d66-8cbf-63a0ed09a2a2

comm_are_nums <- read_csv("raw-data/community_areas_w_numbers.csv", col_names=F) %>%
  select(X1, X2) %>%
  rename(GEOG = X2, area_num_1 = X1) %>%
  mutate(area_num_1 = as.character(area_num_1))
comm_are_nums[comm_are_nums$GEOG=="Lakeview",]$GEOG <- "Lake View"
comm_are_nums[comm_are_nums$GEOG=="Loop",]$GEOG <- "The Loop"

commun_demo <- read_csv("raw-data/CDS_archive_201906/Reference_CCAProfiles_2013_2017.csv") %>%
  select(GEOG, TOT_POP, WHITE, HISP, BLACK, ASIAN, OTHER, MEDINC)

commun_demo <- commun_demo %>%
  left_join(., comm_are_nums, by="GEOG") 
chi <- chi %>% rename(area_num_1=`Community Area`) %>%
  mutate(area_num_1 = as.character(area_num_1))

tot_arrest_data <- chi %>%
  group_by(area_num_1) %>%
  mutate(tot_arrest = sum(Arrest)) %>%
  select(area_num_1, tot_arrest) %>%
  distinct()

tot_arrest_shp <- left_join(commun_area, tot_arrest_data, by="area_num_1") #

tot_arrest_shp$tot_arrest_cat <- cut(tot_arrest_shp$tot_arrest, 
                  breaks = c(0, 500,1000,1500, 2000,2500, 3000, 3500,4000), 
                  labels = c("500","1000","1500", "2000", "2500", "3000", "3500", "4500"))

tot_arrest_shp %>%
  filter(!is.na(tot_arrest_cat)) %>%
ggplot(.) + 
  geom_sf( aes(fill=tot_arrest_cat),size=.3) +
  theme_void() +
  scale_fill_brewer(palette = "RdPu", na.value = "grey80",
                  name="Total Arrests") +
  ggtitle("Chicago Arrests by Community Area, 2018") +
    theme(panel.grid.major = element_line("transparent"), 
        axis.text = element_blank(),
        text=element_text(family="Noto Serif")) 

# join crime data
chi <- chi %>%
  left_join(., commun_demo, by="area_num_1") 
 
chi_shp <- commun_area %>%
  left_join(., chi, by="area_num_1")
# Crime counts by community area
# all reported crimes by community area
tot_crime_data <- chi %>%
  group_by(area_num_1) %>%
  mutate(tot_crime = n()) %>%
  select(area_num_1, tot_crime) %>%
  distinct()

tot_crime_shp <- left_join(commun_area, tot_crime_data, by="area_num_1") #

tot_crime_shp$tot_crime_cat <- with(tot_crime_shp, cut(tot_crime,
   breaks = qu <- quantile(tot_crime, probs = seq(0,1, by=0.2)),
   labels = qu[-1], include.lowest=TRUE))


tot_crime_shp %>%
  filter(!is.na(tot_crime_cat)) %>%
ggplot(.) + 
  geom_sf( aes(fill=tot_crime_cat),size=.3) +
  theme_void() +
  scale_fill_brewer(palette = "YlGnBu", na.value = "grey80",
                  name="Total Crimes") +
  ggtitle("Chicago Crimes Reported by Community Area, 2018") +
    theme(panel.grid.major = element_line("transparent"), 
        axis.text = element_blank(),
        text=element_text(family="Noto Serif")) 

# crime rate per 100,000 people in 2018
crime_rate <- chi_shp %>%
  group_by(area_num_1) %>%
  mutate(crime_rate = (n()/ TOT_POP)*1000) %>%
  ungroup() %>%as.data.frame() %>%
  select(area_num_1, crime_rate) %>%
  distinct()

crime_rate_shp <- commun_area %>%
  left_join(., crime_rate, by="area_num_1") 

crime_rate_shp$crime_rate_cat <- cut(crime_rate_shp$crime_rate, 
                  breaks = c(0, 50,100,150, 200,250, 300, 350), 
                  labels = c("50","100","150", "200", "250", "300", "350"))

ggplot(crime_rate_shp, aes(fill=as.factor(crime_rate_cat))) +
  geom_sf() +
  theme_void() +
    scale_fill_brewer(palette = "YlOrRd", na.value = "grey80",
                  name="Crime Rate") +
  theme(panel.grid.major = element_line("transparent"), 
        axis.text = element_blank(),
        text=element_text(family="Noto Serif"),
        plot.title = element_text(hjust=.3))  +
  ggtitle("Chicago Crime Rate per 1,000 people by Community Area, 2018")

# Black population
comm_dem_shp <- commun_area %>%
  left_join(., commun_demo, by="area_num_1") %>%
  group_by(area_num_1) %>%
  mutate(prop_white = WHITE / TOT_POP,
         prop_black = BLACK / TOT_POP)

ggplot(comm_dem_shp, aes(fill=prop_black)) +
  geom_sf() +
  theme_void() +
  scale_fill_gradient(low = "orange", high = "darkblue", 
                      name = "",
                      labels = percent) +
  ggtitle("Black population by Community Area, 2017") +
    theme(panel.grid.major = element_line("transparent"), 
        axis.text = element_blank(),
        text=element_text(family="Noto Serif"))

311 Light Outages Analysis

library(tidyverse)
# https://data.cityofchicago.org/Service-Requests/311-Service-Requests-Street-Lights-All-Out-Histori/zuxi-7xem/data
allout <- read_csv("raw-data/311_Service_Requests_-_Street_Lights_-_All_Out_-_Historical.csv")

allout$creation_date <-as.Date(allout$`Creation Date`,format="%m/%d/%Y")
allout$completion_date <-as.Date(allout$`Completion Date`,format="%m/%d/%Y")

allout$creation_date_year <- format(allout$creation_date, format = "%Y")
allout <- allout %>% filter(creation_date_year==2018)


# https://data.cityofchicago.org/Service-Requests/311-Service-Requests-Street-Lights-One-Out-Histori/3aav-uy2v/data
oneout <- read_csv("raw-data/311_Service_Requests_-_Street_Lights_-_One_Out_-_Historical.csv")

oneout$creation_date <-as.Date(oneout$`Creation Date`,format="%m/%d/%Y")
oneout$completion_date <-as.Date(oneout$`Completion Date`,format="%m/%d/%Y")
oneout$creation_date_year <- format(oneout$creation_date, format = "%Y")
oneout <- oneout %>% filter(creation_date_year==2018)
lights_out <- allout %>% rbind(.,oneout)
lights_out <- lights_out %>%
  filter(Status=="Completed") %>%
  mutate(time_to_complete_all = completion_date- creation_date) %>%
  group_by(`Type of Service Request`) %>%
  mutate(time_to_complete_status = completion_date- creation_date)
ggplot(lights_out[lights_out$`Type of Service Request`=="Street Lights - All/Out",],
       aes(time_to_complete_status)) + 
  geom_histogram() +
  ggtitle('Distribution of Time to Complete in Days - All/Out') +
  theme_classic() +
  xlab("number of days to complete")

ggplot(lights_out[lights_out$`Type of Service Request`=="Street Light Out",],
       aes(time_to_complete_status)) + 
  geom_histogram() +
  ggtitle('Distribution of Time to Complete in Days - Street Light Out') +
  theme_classic() +
  xlab("number of days to complete")

### Quantiles

# quantile(as.numeric(lights_out[lights_out$`Type of Service Request`=="Street Lights - All/Out",]$time_to_complete_status))
# mean(as.numeric(lights_out[lights_out$`Type of Service Request`=="Street Lights - All/Out",]$time_to_complete_status))
# 
# 
# quantile(as.numeric(lights_out[lights_out$`Type of Service Request`=="Street Light Out",]$time_to_complete_status))
# mean(as.numeric(lights_out[lights_out$`Type of Service Request`=="Street Light Out",]$time_to_complete_status))
# remove na lat/long
lights_out <- lights_out %>%
  filter(!is.na(Latitude))
# convert dataframe to spaital data
lights_out_sf <- st_as_sf(lights_out, coords = c("Longitude", "Latitude"), crs = st_crs(commun_area))

comm_lights_shp <- commun_area %>%
  st_join(lights_out_sf) 

### complete by community area

Light Outages Occurances

comm_lights_shp %>%
  filter(`Type of Service Request`=="Street Lights - All/Out") %>%
  group_by(area_num_1) %>%
  mutate(num_out =n(),
         num_out_cat = cut(num_out,
                 breaks = c(0, 100,200,300,400,500,600), 
                              labels = c("100","200","300","400","500", "600")) 
         ) %>%
  ggplot() +
    geom_sf(aes(fill=as.factor(num_out_cat))) +
    theme_void() +
    scale_fill_brewer(
                  na.value = "grey80",
                  name="Number of All\nLights Out Complaints") +
  ggtitle("Number of 311 All Lights Out Complaints by Community Area, 2018") +
  theme(panel.grid.major = element_line("transparent"), 
        axis.text = element_blank(),
        text=element_text(family="Noto Serif"),
        plot.title = element_text(hjust=.1))

Outages standardized

comm_lights_shp %>%
  left_join(., commun_demo, by="area_num_1") %>%
  filter(`Type of Service Request`=="Street Lights - All/Out") %>%
  group_by(area_num_1) %>%
  mutate(num_out =n(),
         num_out_stand = (num_out/TOT_POP)*1000,
         num_out_stand_cat = cut(num_out_stand,
                 breaks = c(0, 5,10,15,20,25), 
                              labels = c("5","10","15","20","25")) 
         ) %>%
  ggplot() +
    geom_sf(aes(fill=as.factor(num_out_stand_cat))) +
    theme_void() +
    scale_fill_brewer(
                  na.value = "grey80",
                  name="") +
  ggtitle("Number of 311 All Lights Out Complaints per 1,000 by Community Area, 2018") +
  theme(panel.grid.major = element_line("transparent"), 
        axis.text = element_blank(),
        text=element_text(family="Noto Serif"),
        plot.title = element_text(hjust=.5, size=10))

What is the relationship between light outages and crime?

# join light outage data with crime data
crime_lights <- chi %>%
  group_by(Date_simp) %>%
  mutate(crime_n_by_date = n()) %>%
  select(Date_simp, crime_n_by_date) %>%
  rename("creation_date"="Date_simp") %>%
  distinct(creation_date, .keep_all=T) %>%
  right_join(., lights_out, by = c("creation_date"))
##### number of outages on any given day 
lights_out_sub <- lights_out %>%
  ungroup() %>%
  rename(area_num_1 = `Community Area`) %>%
  mutate(area_num_1 = as.character(area_num_1)) %>%
  select(area_num_1, `Service Request Number`, creation_date, completion_date, time_to_complete_status) 

report <- data.frame(date = seq(from = as.Date("2018-01-01"), by="1 day", length.out = 404),
                     count = rep(0,times=404))

c_row <- 0 
for (i in 1:nrow(lights_out_sub)){
  v <- seq(from = as.Date(lights_out_sub$creation_date[i]), by="1 day", 
      length.out = as.numeric(lights_out_sub$time_to_complete_status[i]))
  
  for (j in as.list(v)) {
    report[report$date==j ,]$count <- report[report$date==j ,]$count + 1
  }
  
  c_row <- c_row + 1
  # if (c_row%%1000 == 0){
  #   print(c_row)
  # }
}
##### number of outages on any given day in a community area

final_repot <- data.frame(matrix(ncol = 3, nrow = 0))
for (area_num in unique(lights_out_sub$area_num_1)) {
  # print(area_num)
  
  lo_grp <- lights_out_sub %>%
    filter(area_num_1==area_num)
  
  report_comm <- data.frame(date = seq(from = as.Date("2018-01-01"), by="1 day", length.out = 404),
                     count = rep(0,times=404),
                     area_num = rep(area_num,times=404))
  
  for (i in 1:nrow(lo_grp)){
    v <- seq(from = as.Date(lo_grp$creation_date[i]), by="1 day", 
        length.out = as.numeric(lo_grp$time_to_complete_status[i]))

    
        for (j in as.list(v)) {
          report_comm[report_comm$date==j ,]$count <- report_comm[report_comm$date==j ,]$count + 1
        }
        
        # rbind report at the end of group loop
        if (i==max(nrow(lo_grp))) {
          final_repot <- rbind(final_repot,report_comm)
        }

  }
}


# st_write(commun_area, "./output-data/commun_area.shp")
# average number of light outages in July by community area
july_val <- final_repot %>%
  mutate(date = as.Date(date),
         month = month(date)) %>%
  filter(month==7) %>%
  group_by(area_num) %>%
  summarise(m = mean(count)) %>%
  rename(area_num_1 = area_num) 
# average number of light outages in July by community area
july_val <- final_repot %>%
  mutate(date = as.Date(date),
         month = month(date)) %>%
  filter(month==7) %>%
  group_by(area_num) %>%
  summarise(m = mean(count)) %>%
  rename(area_num_1 = area_num) 

# number of light outages in July by community area --- standardized
july_val_pop <- july_val %>%
  left_join(., commun_demo, by="area_num_1") %>%
    mutate(num_out =n(),
         num_out_stand = (num_out/TOT_POP)*1000,
         num_out_stand_cat = cut(num_out_stand,
                 breaks = c(0, 5,10,15,20,25), 
                              labels = c("5","10","15","20","25")) 
         ) 

july_val_pop$num_out_stand_cat_quant <- with(july_val_pop, cut(num_out_stand,
   breaks = qu <- quantile(num_out_stand, probs = seq(0,1, by=0.2), na.rm=T),
   labels = round(qu[-1],2), include.lowest=TRUE))

commun_area %>%
  left_join(., july_val_pop, by="area_num_1") %>%
  ggplot(., aes(fill=as.factor(num_out_stand_cat_quant))) +
  geom_sf() +
  theme_void() +
    scale_fill_brewer(
                  na.value = "grey80",
                  name="") +
  ggtitle("Number of Incomplete All Lights Out Complaints per 1,000 by Community Area, July 2018") +
  theme(panel.grid.major = element_line("transparent"), 
        axis.text = element_blank(),
        text=element_text(family="Noto Serif"),
        plot.title = element_text(hjust=.5, size=10)) 

# crime rate in July by community area --- standardized

crime_rate_july <- chi_shp %>%
  filter(month==7) %>%
  group_by(area_num_1) %>%
  mutate(crime_rate = (n()/ TOT_POP)*1000) %>%
  ungroup() %>%as.data.frame() %>%
  select(area_num_1, crime_rate) %>%
  distinct()

crime_rate_shp_july <- commun_area %>%
  left_join(., crime_rate_july, by="area_num_1") 

crime_rate_shp_july$crime_rate_cat_quant <- with(crime_rate_shp_july, 
                                                    cut(crime_rate,
   breaks = qu <- quantile(crime_rate, probs = seq(0,1, by=0.2), na.rm=T),
   labels = round(qu[-1],2), include.lowest=TRUE))

ggplot(crime_rate_shp_july, aes(fill=as.factor(crime_rate_cat_quant))) +
  geom_sf() +
  theme_void() +
    scale_fill_brewer(palette = "YlOrRd", na.value = "grey80",
                  name="Crime Rate") +
  theme(panel.grid.major = element_line("transparent"), 
        axis.text = element_blank(),
        text=element_text(family="Noto Serif"),
        plot.title = element_text(hjust=.3))  +
  ggtitle("Chicago Crime Rate per 1,000 people by Community Area, July 2018")

# average night crime rate in July by community area --- standardized

crime_rate_july <- chi_shp %>%
  mutate(time = hms::as_hms(time)) %>%
  filter(month==7,
         time >=hms::as_hms('21:00:00') | time <=hms::as_hms('6:00:00')
         ) %>%
  group_by(area_num_1) %>%
  mutate(crime_rate = (n()/ TOT_POP)*1000) %>%
  ungroup() %>%as.data.frame() %>%
  select(area_num_1, crime_rate) %>%
  distinct()

crime_rate_shp_july <- commun_area %>%
  left_join(., crime_rate_july, by="area_num_1") 

crime_rate_shp_july$crime_rate_cat_quant <- with(crime_rate_shp_july, 
                                                    cut(crime_rate,
   breaks = qu <- quantile(crime_rate, probs = seq(0,1, by=0.2), na.rm=T),
   labels = round(qu[-1],2), include.lowest=TRUE))

ggplot(crime_rate_shp_july, aes(fill=as.factor(crime_rate_cat_quant))) +
  geom_sf() +
  theme_void() +
    scale_fill_brewer(palette = "YlOrRd", na.value = "grey80",
                  name="Crime Rate") +
  theme(panel.grid.major = element_line("transparent"), 
        axis.text = element_blank(),
        text=element_text(family="Noto Serif"),
        plot.title = element_text(hjust=.3))  +
  ggtitle("Chicago Nightly Crime Rate per 1,000 people by Community Area, July 2018")

## daily nightly crime rate
## daily light outages per 1,000

crime_rate_july_daily <- chi_shp %>%
  mutate(time = hms::as_hms(time)) %>%
  filter(month==7,
         time >=hms::as_hms('21:00:00') | time <=hms::as_hms('6:00:00')
         ) %>%
  group_by(area_num_1,Date_simp) %>%
  mutate(crime_rate = (n()/ TOT_POP)*1000) %>%
  ungroup() %>%as.data.frame() %>%
  select(area_num_1, Date_simp, crime_rate) %>%
  distinct()

daily_july_out_stand <- final_repot %>%
  rename(area_num_1 = area_num) %>%
  left_join(., commun_demo, by="area_num_1") %>%
  group_by(area_num_1,date) %>%
    mutate(num_out_stand = (count/TOT_POP)*1000,
           date = as.Date(date),
         month = month(date)) %>%
  filter(month==7) %>%
  select(area_num_1, date, num_out_stand)

the_data <- crime_rate_july_daily %>%
  rename(date = Date_simp) %>%
  left_join(., daily_july_out_stand, by=c("area_num_1", "date"))

GWR: Geographically Weighted Regression

library(spgwr)

gwr_df <- commun_area %>%
  left_join(., the_data, by="area_num_1")

gwr_df_cent <- st_centroid(gwr_df)
## Warning in st_centroid.sf(gwr_df): st_centroid assumes attributes are constant
## over geometries of x
## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for longitude/
## latitude data
gwr_df_cent <- gwr_df_cent %>%
  mutate(seq_day = day(date)
)

GWRbandwidth <- gwr.sel(crime_rate~num_out_stand+as.factor(seq_day),
                        data=gwr_df_cent, 
                        coords=st_coordinates(gwr_df_cent),adapt=T)
## Adaptive q: 0.381966 CV score: 19.25458 
## Adaptive q: 0.618034 CV score: 20.12083 
## Adaptive q: 0.236068 CV score: 18.37588 
## Adaptive q: 0.145898 CV score: 17.5352 
## Adaptive q: 0.09016994 CV score: 16.93422 
## Adaptive q: 0.05572809 CV score: 16.36893 
## Adaptive q: 0.03444185 CV score: 15.85204 
## Adaptive q: 0.02128624 CV score: 15.94012 
## Adaptive q: 0.03158602 CV score: 15.91558 
## Adaptive q: 0.04257247 CV score: 16.03327 
## Adaptive q: 0.03754747 CV score: 15.91985 
## Adaptive q: 0.03451831 CV score: 15.85166 
## Adaptive q: 0.03490357 CV score: 15.8406 
## Adaptive q: 0.03591345 CV score: 15.81673 
## Adaptive q: 0.03653759 CV score: 15.93889 
## Adaptive q: 0.03552771 CV score: 15.83346 
## Adaptive q: 0.03615185 CV score: 15.84553 
## Adaptive q: 0.03576611 CV score: 15.82167 
## Adaptive q: 0.03600451 CV score: 15.8152 
## Adaptive q: 0.0360452 CV score: 15.82251 
## Adaptive q: 0.03596382 CV score: 15.81535 
## Adaptive q: 0.03600451 CV score: 15.8152
#run the gwr model
gwr.model = gwr(crime_rate~num_out_stand+as.factor(seq_day),
                        data=gwr_df_cent, 
                        coords=st_coordinates(gwr_df_cent),
                adapt=GWRbandwidth, hatmatrix=TRUE, se.fit=TRUE) 

#extract results for each variable
results<-as.data.frame(gwr.model$SDF)
#attach coefficients to original dataframe
gwr_df_cent$coefest_num_out_stand<-results$num_out_stand

gwr_df_cent$coefbls_as.factor.seq_day.10<-results$as.factor.seq_day.10

summary(results$num_out_stand) #same summary data as you see in gwr.model object
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.08195  0.02094  0.09737  0.11204  0.20463  0.35390
#now plot the various GWR coefficients

# What is the relationship between community area daily light outages and daily crime rate across Chicago? Is the relationship uniform?
library(leaflet)
library(leaflet.providers)

qpal<-colorQuantile("OrRd", gwr_df_cent$coefest_num_out_stand, n=9) 

st_crs(commun_area) <- '+proj=longlat +datum=WGS84'
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that
st_crs(gwr_df_cent) <- '+proj=longlat +datum=WGS84'
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that
leaflet(gwr_df_cent) %>%
  addTiles() %>%  
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(data=commun_area, color = "black", fillOpacity = 0, opacity=1, weight=1) %>%
  addCircleMarkers(radius=3,color = ~qpal(gwr_df_cent$coefest_num_out_stand)
  ) %>%
  addLegend(position = "bottomright",
              pal = qpal, values = ~gwr_df_cent$coefest_num_out_stand,
              title = "Standardized Lights Out Residuals")

The relationship between daily number of lights out per 1,000 and crime rate is stronger, or has a larger effect size, in the central-sourthern region of the city